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Abstract 


Thermal-fluid transport phenomena in a membrane-electrode assembly (MEA) of a polymer electrolyte membrane (PEM) fuel cell attached to 
interdigitated gas distributors are studied numerically. The MEA consists of two porous catalyst layers, two porous gas diffusion layers, and an 
impermeable PEM. In the catalyst layers, the overpotential heating by the electrochemical reaction under thermal equilibrium conditions produces 
heat that is removed by the fluids as well as the solid matrices. In the diffusion layers, the difference in the heat conductivities between the solid 
matrices and the fluids causes a thermal non-equilibrium in the porous medium. A two-equation approach is used to resolve the temperature difference 
between the solid matrices and the fluids. The effects of the porous Reynolds number, interfacial heat transfer coefficient, and overpotential heating 
are examined. It is found that the local maximum temperature occurs inside the cathodic catalyst layer. In addition, the temperature difference 
between the solid matrices and the fluids in the diffusion layers decreases with increasing the non-dimensional interfacial heat transfer coefficient. 
The present results have provided comprehensive heat transfer information that is helpful in understanding of the mechanisms responsible for 


thermal pathways in a PEM fuel cell. 
© 2006 Elsevier B.V. All rights reserved. 
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1. Introduction 


Fuel cells convert directly the chemical energy of fuels 
to electricity through electrochemical reactions. The combi- 
nation of benefits in high efficiency, eco-friendly nature, and 
diversity in fuels makes them attractive alternative power con- 
version devices. The classification of fuel cells is commonly 
based on the electrolyte material, including polymer electrolyte 
membrane (PEMFC), solid oxide (SOFC), molten carbonate 
(MCFC), potassium hydroxide (AFC), and phosphoric acid 
(PAFC). Among them, PEM fuel cells have appeared as the 
most promising candidates for commercial applications. 

To improve the understanding of the characteristics in a PEM 
fuel cell, a large amount of experimental efforts have been 
devoted in the past decade. Common techniques include the 
ac-impedance approach that determines the impedance char- 
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acteristics and accesses the modes of failure of a PEM fuel 
cell, and the polarization-curve measurement that provides the 
operating J-V data of a PEM fuel cell. In general, these experi- 
ments provide overall information but lack for revealing detailed 
characteristics inside a PEM fuel cell such as local temperature 
distributions. It is known that the efficiency and durability of a 
PEM fuel cell are strongly affected by the temperature distribu- 
tions inside the fuel cell. The present instruments cannot support 
such measurements due to the highly reactive environment and 
compactness of a fuel cell. Thus, thermal-fluid transport phe- 
nomena in the PEM fuel cell are still not well understood. Such 
information should be sought relying on modeling or simulation. 

A large amount of fuel cell models have been developed 
accounting for various physical processes. Due to space lim- 
itation, only some relevant works are briefly reviewed below. 
Bernardi and Verbrugge [1] proposed a one-dimensional com- 
putational model to address the water management and species 
transport in the cathodic gas channel and gas diffusion layer 
attached to the membrane. Assumptions used in this model are 
isothermal, steady state, and ideal gases. Springer et al. [2] pre- 
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Nomenclature 


Bi Biot number, Table 2 

CHO water vapor mole concentration (mol m~?) 

CH)O, ref Water vapor mole concentration at inlet 
(mol m~?) 

Cmo normalized water vapor concentration, Table 2 

Cy50, ref normalized water vapor mole concentration at 
inlet, Table 2 


co, oxygen concentration (mol m7?) 

CO>,ref Oxygen mole concentration at inlet (mol m7?) 

Co, normalized oxygen concentration, Table 2 

Co,ref normalized oxygen mole concentration at inlet, 
Table 2 

Cp specific heat at constant temperature (J kg~! K7!) 

Ctot total mole concentration of the reacting fluid 
(mol m~?) 

d pore diameter of the porous medium (m) 


Dp, hydrogen diffusivity (m? s7!) 

Dy, ett effective diffusivity of hydrogen (m? s7!) 
Do, diffusivity of oxygen (m? s7!) 

Do,,eff effective diffusivity of oxygen (m? s7!) 
Da Darcy number, Table 2 


Eco,, Ec, , Ecm o non-dimensional electrochemical- 
reaction heating, Table 2 

F Faraday’s constant (96487 C mol— 9) 

hy interstitial heat transfer coefficient (W m~? K7!) 

Jo exchange current density (A m7?) 

JT transfer current density (A m~>) 

k thermal conductivity (W m7! KE 


Mto,, Mty,, Mty,0 non-dimensional mass transfer coef- 
ficient, Table 2 


p pressure (Pa) 

P non-dimensional pressure, Table 2 

Pr Prandtl number, Table 2 

Rk ratio of fluid-to-solid conductivity, Table 2 
Re Reynolds number, Table 2 

Se source term for the energy equation 

Si source term for the concentration equation 
Sm source term for the momentum equation 
Sc Schmidt number, Table 2 

T temperature (K) 

u, U velocity vectors 

uq pore velocity in the porous medium (m s7!) 


U,V non-dimensional velocity components in the X, 
and Y direction 

W witdth of the computational domain (m) 

X, Y coordinate system, Fig. 2 


Greek letters 
œ1,œ2 coefficients in Eq. (9) 


Bi coefficient in Eq. (8) 

ô thickness of the five-layer MEA (m) 
E porosity 

n overpotential (V) 

K permeability (m7) 


0 non-dimensional temperature, Table 2 
p density (kg m~?) 

T tortuosity of the porous electrode 
Subscripts 

a anode 

c cathode or catalyst layer 

CL catalyst layer 

d diffusion layer 

eff effective 

f fluid phase 

HOR hydrogen oxidation reaction 

i species 

m membrane 

ORR oxygen reduction reaction 

ref reference 

s solid phase 

tot total 


sented an isothermal, one-dimensional, steady state model for 
a PEM fuel cell with a hydrated Nafion-117 membrane. The 
unique feature of this model is the empirical relation for cal- 
culating membrane conductivity based on water content in the 
membrane. Nguyen and White [3] presented a steady state, two- 
dimensional mass transfer model. The model domain consists of 
gas flow channels, gas diffusion layers, catalyst layers, and mem- 
brane. Gurau et al. [4] applied computational fluid dynamics to 
fuel cell modeling where they began with a two-dimensional 
model consisting of both the anode and cathode sides. This 
model fully accounted for the mass transport of reactant species 
which have been simplified in earlier models. Yi and Nguyen [5] 
modified their previous model to include heat and mass trans- 
fer conditions in both the liquid and gas phases along the flow 
path of the fuel cell. Shimpalee and Dutta [6] described a steady 
state, isothermal, three-dimensional, and single phase PEM fuel 
cell model. A commercial code was used to solve the complete 
Navier-Stokes equations. Hwang et al. [7—9] described a non- 
isothermal model for the cathodic half-cell of a PEM fuel cell. 
Non-dimensional forms of the governing equations were devel- 
oped and the effects of dimensionless parameters on fuel cell 
performance were evaluated. 

From the above discussion, it is obvious that most previ- 
ous fuel cell models regarded the electrochemical reaction as 
an isothermal process. Studies of thermal-fluid transports in the 
PEM fuel cells are rather sparse. The present study extends the 
authors’ previous work of cathodic heat transfer [7—9] to study 
the heat/mass transfer in an entire membrane-electrode assem- 
bly of a PEM fuel cell. A computational model considering local 
thermal non-equilibrium (LTNE) in the gas diffusion layers are 
developed to study the detailed thermal-fluid behaviors in a PEM 
fuel cell. Although, the LTNE analysis of heat transfer in the 
porous medium is not new [10-14], a systematical study on the 
LTNE characteristics of a five-layer membrane-electrode assem- 
bly of a PEM fuel cell that is coupled with the electrochemical 
reaction might be considered to be original and valuable. 
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2. Numerical model 
2.1. Model descriptions 


A typical PEM fuel cell with interdigitated gas distributors 
is schematically shown in Fig. 1. The computational domain 
is confined within a five-layer membrane-electrode assembly 
(MEA), i.e., two gas diffusion layers (GDL), two catalyst layers 
(CL), and a polymer electrolyte membrane (PEM). Electro- 
chemical reactions in the MEA include the hydrogen oxida- 
tion reaction (HOR) in the anodic catalyst layer and the oxy- 
gen reduction reaction (ORR) in the cathodic catalyst layer. 
The hydrogen enters the anode from the bottom-left corner 
of the module, traverses the anodic diffusion layer to the 
anodic catalyst layer, and is consumed by HOR. Meanwhile, 
the protons and heat are produced in the anodic catalyst layer, 
i.e., 


5H: > H+ + e7 + Heat (1) 


The un-reacted hydrogen exits by the bottom-right corner of the 
module. In the cathode, the oxygen enters the cathodic diffusion 
layer from the upper-left corner of the module, transverses the 
diffuse layer to the catalyst layer, and then reacts together with 
the protons to form water and heat. 


102 +H* +e7 + 5H20 + Heat (2) 


Finally, the products (water vapor) exit by the upper-right corner 
of the module. 

Physical and geometric parameters of each layer of the MEA 
are given in Table 1. The coordinate system of the computational 


Ozin O2, H20 out 


Cathodic Catalyst Layer 
(Porous Medium) 


Anodic Catalyst Layer 
(Porous Medium) 


Hout 


Computational Domain 
5-layer Membrane-Electrode Assembly 


Oxygen Inlet Channel Oxygen Outlet Channel 


Interdigitated 
Gas Distributor 


Hydrogen Inlet Channel Hydrogen Outlet Channel 


Fig. 1. Schematic drawing of the five-layer MEA of a PEM fuel cell with inter- 
digitated gas distributors. 


Table 1 
Physical and geometric parameters of the five-layer MEA 
Description Value 
Porosity, € 
Gas diffusion layer, EGDL 0.48 
Catalyst layer, ec. 0.42 


Permeability, « (m?) 


Gas diffusion layer, KGDL 1.57 x 107!2 
Catalyst layer, kcL 1.02 x 107!2 
Tortuosity, T 
Gas diffusion layer, TGDL 1.5 
Catalyst layer, TCL 1.5 
Thermal conductivity, k (Wkg—! K7!) 
Gas diffusion layer, ks eff 1,7 
PEM, km 0.5 
Anode gas, kf ert 0.182 
Cathode gas, kf,eff 0.051 
MEA dimensions (pm) 
Gas diffusion layer thickness, dgpL 200 
Catalyst layer thickness, dcL 25 
PEM thickness, ôm 50 
MEA thickness, ô 500 


domain is shown in Fig. 2. The assumptions used in this model 
are as follows [15]: 


(1) Gas mixtures are ideal gas. 

(2) The fluid flow is steady, laminar, and incompressible; its 
thermal physical properties are constant. 

(3) Porous electrode is homogeneous and isotropic with uni- 
form morphological properties such as porosity, tortuosity 
and permeability. 

(4) Water in the electrode exits as vapor only. 

(5) Catalyst layer is treated as an ultra-thin layer; thus the oxy- 
gen reduction reaction is considered to occur only at the 
surfaces of the catalyst layer. 

(6) The inlet fluid and rib-surface temperatures are uniform. 


0.4 0.1 1.6 


Fig. 2. Dimensions and coordinate system of the computational module. 
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2.2. Momentum transports 
The mass conservation equation together with the Brinkman- 


extended Darcy equation describes the fluid flow in the gas 
diffusion layers. 


V(pu) = 0 (3) 


pu: Vu = —Vp + V : (Vu) + Sm (4) 


The source term in the momentum equations is Sm = —(u/K)eu, 
respectively. Using the non-dimensional parameters list in 
Table 2, the above equations are further reduced to dimension- 
less forms of: 


VU=0 (5) 

U-vu=-ve+ vu ê U] (6) 
~ Re Re- Da 

Table 2 


Non-dimensional parameters used in the present study 


2.3. Concentration distributions 


The species transports in the electrodes are described by the 
following equations: 


u: Vc = Di V’ ci + Si (7) 


The species considered in the anode and cathode are pure hydro- 
gen, cp, and the wet oxygen, co, and cp,o, respectively. The 
effective diffusivities of the species i in the porous electrode 
follows the Bruggemann model [16], i.e., Di eff = e'Di. The 
source term is nothing in the gas diffusion layers. In the cat- 
alyst layer, it is implemented based on electrochemical kinetics, 
i.e., consuming the reactants to produce the current. In the 
anodic catalyst layer, it becomes $;=—(tHoR/2F). Accord- 
ing to the Butler—-Volmer correlation [17], the rate of electro- 
chemical reaction on the reaction surfaces can be described by 
the relationship of the local current density and the reactant 


Symbols Expression Physical meanings 
: hy y 
Bi — Biot number 
ks eff 
CH, Cirat i $ P 
Cup, Cup, ref a ee Non-dimensional hydrogen concentration 
Ctot  Ctot 
CO, €Op, ref : ; p 
Co, Coy, ref a2 Non-dimensional oxygen concentration 
Ctot  Ctot 
K 
Da z Darcy number 
-8n f , , ee : 
Ecy, Pi li Ratio of the electrochemical heating to the conduction in solid phase of the anode 
ks,eff(Tf,in a T:)Ch; ref 
q2: & “Ne ; P : ‘ ‘ : . 
Ecm, o z Ratio of the electrochemical heating to the conduction in solid phase of the cathode 
ks eff(Tf,in ai Te)Cho,ref 
ay: 82. Ne : P Í Ki a z 
Eco, Ratio of the electrochemical heating to the conduction in solid phase of the cathode 
ks eff(Tf,in = T;)COp ref 
Joa: Bi 8 ' ' PERE 
Mty, — mM Ratio of the hydrogen consumption to the mass diffusion in the anode 
. 2F + CHy ret * Dib eff 
Jo c' a2: & A è 7 $ F 
Mty,0 E Ratio of the hydrogen consumption to the mass diffusion in the anode 
4F + ciot: CHO, ref Dy) eff 
joc: a18? 
Mto, eS Ratio of the oxygen consumption to the mass diffusion in the anode 
a 4F - Coy ref © Doo eft 
P za Non-dimensional pressure 
prs 
(Cp) 
Prese Pat Prandtl number 
kr eff 
ugd 
Re — Reynolds number 
v 
Rk feti Ratio of the fluid-phase conductivity to the solid-phase conductivity 
ks eff 
v v 
Scy,, Sco, =; Schmidt numbers for hydrogen and oxygen, respectively 
Duz et Donet 
u v : ; : 
U,V U=—,V= — Non-dimensional velocity 
ud ud 
T-T : A 
6 : Non-dimensional temperature 


Ty.in al Te 


454 J.J. Hwang et al. / Journal of Power Sources 163 (2006) 450-459 


concentrations, 1.e., 


2 . CH 
JT,HOR = Jjo,aPi ( 2 ) (8) 
CH), ref 
co cmo \? 
JT,OOR = Jo,c |&1 ( 2 ) a( 2 ) (9) 
CO> ref CH) O,ref 


where jo,a and jo, are the anodic and cathodic exchange cur- 
rent densities, respectively. a1, a2, and f; are electrochemical 
coefficients depending on and the overpotential on the electrode 
surfaces. All above coefficients are regarded as constants in the 
present simulation. Based on the non-dimensional parameters in 
Table 2, the dimensionless forms of the species transport equa- 
tions in the anodic and cathodic catalyst layers are, respectively, 
written as 


2 M 2 
U-VCu, = V?Cm + —— Cp, (10) 
R H 


U-VCo, = 


| 
N 
Q 

© 


— [Mts Co, = Mtp-o(l — Co,)] (1 
Re Seo," o Co, H,0( oy] dt) 


2.4. Temperature distributions 


In practical applications, the temperature difference between 
the inlet and the outlet of a PEM fuel cell is not so large, typi- 
cally about 30-40 K. In addition, the conductivities of the solid 
matrix (such as carbon fibers) and the reactant fluid (such as air) 
are quite different. That is the temperature between the solid 
matrices and the fluids may be different and thus away from 
the local thermal equilibrium. Therefore, in the present model, 
a two-equation model is used to describe the thermal behaviors 
in the gas diffusion layers, i.e., 


(pcp): VTe = V - (kt eteV Te) — Se,GDL (12) 
0=V. (ks, ett V Ts) F Se, GDL (13) 


The source terms in the above equations represent the thermal 
interaction between the solid matrices and the fluids. They are 
represented by Se Got = —/y - (T; — Te), where hy is the inter- 
facial heat transfer coefficient (volumetric) between the solid 
matrices and the reactant fluids in porous medium. Actually, up 
the present time, the values of hs for porous electrode has not 
been determined yet, and thus is not available in the open litera- 
ture. According to the measured data of the metal foams [11], the 
typical values of hy vary from 2.0 x 104 to 2.0 x 105 W m~? K7! 
for the porosity varying from 0.7 to 0.95. In the present study, 
the value of hy for the baseline case (St= 1.39) is set to be 
1.0 x 10° W m~? K- !. The effective thermal conductivities of the 
solid phase and the fluid phase are, respectively, defined as 


ks eff = (1 — £)ks (14) 


kf eff = ekg (15) 


where ks and kg are the thermal conductivities of the carbon fiber 
and the reactant fluid, respectively. 

In the catalyst layers, the electrochemical reaction for the 
reactant fluids should take place on the catalyst surfaces under 
the same temperature (i.e., reacting temperature). Thus, the flu- 
ids and the solid matrices should be in thermal equilibrium. 


T: = Ts (16) 
(pcp) U- VT; = V - (kcL VT) + Sect (17) 


In the present model, the radiation heat flux and the energy dissi- 
pation due to Joule heating are neglected. Therefore, the source 
term in the catalyst layer is represented by the overpotential heat- 
ing by the electrochemical activation, i.e., Se CL =JTHOR X Ma in 
the anodic catalyst layer and Se,CL =JT,ORR X Nc in the cathodic 
catalyst layer, respectively. 

The effective thermal conductivity of the catalyst layer is 
determined by the following equation: 
2ke + : (18) 

© (€/(Qke + ks) + (1 = €)/3ke) 
where k, is the weight-averaged conductivity between the ionic 
conductor (such as Nafion™) and the electric conductor (such 
as Pt/C). 

Using the non-dimensional parameters given in Table 2, the 
dimensionless forms of the energy equations in the diffusion 
layers for both the anode and the cathode are 


kcL = 


1 2 
U- V6 = ——_V°6 + 


Bi 
ih pe 19 
Re. Pr Re Pro R, n a 


0 = V70; — Bi(Os — Of) (20) 


where Bi is the non-dimensional interfacial heat transfer coeffi- 
cient (Table 2). In general, it is a function of morphology of the 
porous matrices. Increasing the surface-to-volume ratio of the 
porous matrices can enhance hy that reduced the temperature 
difference between the solid matrix and the reactant fluid. 

In the cathodic and anodic catalyst layers, the dimensionless 
energy equations become 


U. Ve = V7O¢ 
Re. Pr 
2 
tgz. py E0200 — Eemol — Co,)"] (21) 
U- V6 = -Vet tc (22) 
= Re- Pr © Re Pr 


In addition, the fluid and the solid have the same temperatures, 
i.e., 


s = Of (23) 


As for the PEM, a typical conduction equation is employed 
to describe the thermal behavior in the impermeable material, 
L.e., 


V - (km VT.) = 0 (24) 
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0.0020 
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0.0004 £7" 
0.0002 Ë 
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Fig. 3. Effect of Reynolds number on the flow velocity distribution in the MEA. 


where km is the ionic conductivity of the PEM. The non- 
dimensional form of the conduction equation in the PEM 
becomes 


ver =0 (25) 


2.5. Boundary conditions 


The temperatures for both feeds along with the surfaces of the 
gas distributor are constant, i.e., 0 = 1. Both outlets of the module 
have a reference (ambient) pressure. The anode is supplied with 
pure hydrogen, while the cathodic side feeds with the humidified 
oxygen of 0.9/0.1 for Co, /Cm,o. Flow rates of the reactants 
delivered to the cathode and anode are varied by examined the 
Reynolds number effect (Re = 1-20). 


2.6. Numerical methods 


In the present study, the finite element solutions are obtained 
with a general purpose commercial solver, COMSOL Multi- 
1.00 
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Fig. 5. Distributions of local minimum concentration of hydrogen and oxygen 
along the flow direction. 


physics [18-22]. It uses the Broyden’s method with an LU- 
decomposition pre-conditioner to solve the non-linear equa- 
tions iteratively. To reduce continuity errors, a penalty term 
is employed for pressure. Thus, there is a continuous part of 
the pressure and piecewise constant part providing and extra 
DOF (degree of freedom) for pressure on each element. It uses 
Newton-Raphson iteration to solve the close-coupled groups 
(velocity, pressure, temperature, concentration and electricity) 
and uses the frontal algorithm (Gaussian elimination) to solve 
the linearized system of equations for each iteration. 

In the present work, all computations are performed on 
50 x 160(X x Y) structured, orthogonal meshes. Additional runs 
for the coarser meshes, 40 x 120, and the finer meshes, 60 x 180, 
are taken for a check of grid independence. A comparison of the 
results of the two grid sizes, 50 x 160, and 60 x 180, shows that 
the maximum discrepancies in the axial velocity and oxygen 
concentration profiles are only 0.6 and 0.9%, respectively. In 
addition, results indicate a maximum change of 0.6% in wall 
temperature distribution between the solutions of 50 x 160 and 


1.00 
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0.90 
0.88 2 
0.86 
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Fig. 4. Species concentration distributions at several axial stations (X) in the MEA. 


456 


60 x 180 grids. All above discrepancies are so small that the 
accuracy of the solutions on a 50 x 160 grids is deemed satis- 
factory. A typical simulation requires about 200 min of central 
processing unit time on a Pentium IV 2.8 GHz PC. 


3. Results and discussion 


Before the discussion of the numerical results, it requires to 
validate the numerical model by comparing the present numeri- 
cal results with the available experimental data. Up to the present 
time, it is hard to measure the local temperature inside the porous 
electrode of a fuel cell. Strictly speaking, no data about the 
local thermal-fluid characteristics such as local temperatures 
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and heat transfer coefficients inside the porous electrode of a 
fuel cell are available in the open literature. Therefore, a macro- 
scopic comparison of polarization curve is made between the 
present predictions and previous experiments [23], which have 
been widely used to verify a fuel cell simulation model. A com- 
parison of the polarization curve between the present numerical 
predictions and the experimental data under adiabatic conditions 
showed a satisfactory agreement. 

Fig. 3 shows the flow velocity distribution (vV U? + V?) at 
the module middle (X = 0.8) as a function of the Reynolds num- 
ber. The Reynolds number varies from Re = 1 to 20. The Darcy 
numbers are different between the catalyst layer and the gas 
diffusion layer due to the difference in the permeability. They 
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Fig. 6. Comparison of fluid-phase and solid-phase temperature distributions at several X stations in the MEA for (a) Bi/(Pr-Re-R,)=6.63 x 107°, (b) 


BiKPr-Re-Ry) =0.663, (c) Bi/(Pr-Re-Ry) = 6.63, and (d) Bi/(Pr-Re-R;) = 6.63 x 10°. 
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are Da=4.9 x 1076 and Da=9.8 x 107° for the catalyst layers 
and the gas diffusion layers, respectively. No data are shown in 
the membrane (0.45 < Y < 0.55) since it is impermeable. The gas 
diffusion layers have a higher permeability and a higher porosity 
higher than the catalyst layers. Therefore, a higher velocity is 
accompanied by the gas diffusion layers. In addition, the maxi- 
mum velocities occur near the shoulder surfaces (Y=0 and 1.0) 
for the both electrodes. It is reasonable because the pathways 
from the module inlet to the module outlet are shorter for the 
flow closer to the shoulder surfaces. A shorter flow pathway pays 
a smaller pressure-drop penalty, and thus allows more fluids to 
pass through. It is further seen that the velocity increases with 
increasing the Reynolds number in both the catalyst layers and 
the diffusion layers. In addition, the flow velocity distributions 
show small humps near the interfaces of the catalyst layer and 
the gas diffusion layer (Y= 0.4 and 0.6) and near the rib shoulder 
(Y=0 and 1.0) due to the shear flows. 

Fig. 4 shows the concentration distributions of oxygen (Co, ) 
and hydrogen (Cy, ) at several axial (X) stations of the porous 
electrodes. The species concentrations are Cy, = 1.0 (pure 
hydrogen) for the anode inlet and Co, = 0.9 and Cp,o = 0.1 
(humidified oxygen) for the cathode, respectively. The Reynolds 
number is fixed at Re = 6. At the axial station of X = 0.4, the pure 
hydrogen decreases its concentration from the inlet to the cat- 
alyst layer. As the flow moves downstream, Cy, is gradually 
decreased due to the HOR in the catalyst layer. A similar trend 
is also observed for the oxygen concentration in the cathode. The 
lowest values of Cy, and Co, at each axial station occur at the 
interfaces of the membrane and the catalyst layers, i.e., Y= 0.45 
and 0.55, respectively. Fig. 5 further shows the distributions of 
the local minima of Cy, and Co, along the axial direction (X). 
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As the flow goes downstream, the distributions of the local min- 
ima of Cy, and Co, start with a little decline, then drop sharply, 
and finally decrease mildly. 

Fig. 6 compares the distributions of the fluid-phase tempera- 
ture (0f) and the solid-phase temperature (6,) in the MEA under 
various Bi/(Pr-Re-R;). Again, no fluid-phase temperatures are 
shown in the membrane due to its impermeability. In the cata- 
lyst layers (0.4 < Y < 0.45, 0.55 < Y<0.6), the solid matrices and 
the fluids have the same temperatures. A local maximum tem- 
perature is clearly observed at about Y= 0.57, where significant 
heat is generated by the electrochemical reaction. In the diffusion 
layer, the fluid-phase temperature and solid-phase temperature 
separate from the interface of the catalyst layer and the gas dif- 
fusion layer due to the LTNE effect. The extent of temperature 
discrepancy between these two phases in the gas diffusion layer 
depends on Bi/(Pr-Re-R,) as well as the location. At the smallest 
value of Bi/(Pr-Re-R;,) = 6.63 x 107? (Fig. 6(a)), the fluid-phase 
temperature and the solid-phase temperature are rather different 
at the axial stations of X=0.8 and 1.2. It is further seen that 
an increase in the value of Bi/(Pr-Re-R;) reduces the tempera- 
ture discrepancy between the solid phase and the fluid phase. As 
the value of Bi/(Pr-Re-R;,) increases to 6.63 x 107? (Fig. 6(d)), 
the fluid-phase temperature and the solid-phase temperature fall 
into a single curve almost. A great similarity in the temperature 
distribution between the solid phase and the fluid phase gives 
an indication that the thermal-fluid field in the porous electrode 
has become equilibrium locally. 

Fig. 7(a and b) compare the temperature distributions in the 
MEA between two different values of Bi/(Pr-Re-R;) at three 
axial stations of X = 0.4, 0.8 and 1.2. As shown in Fig. 7(a) 
for the smaller value of Bi/(Pr-Re-R,;) (6.63 x 10~3), the fluid- 
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Fig. 7. Fluid-phase and solid-phase temperature distributions at several X stations in the MEA for (a) Bi/(Pr-Re-Rx) = 6.63 x 1073, and (b) Bi/(Pr-Re-R,) = 6.63. 
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Fig. 8. Temperature distributions along X direction at several Y stations for Bi/(Pr-Re-Rx) = 6.63, (a) solid-phase temperature, and (b) fluid-phase temperature. 


phase temperature increases monotonically downstream in both 
electrodes. The solid-phase temperature decreases downstream 
in the core region of the MEA but has the lowest value near 
the middle of the shoulder surfaces (X = 0.8). This is because 
the module middle has the shortest distance for heat conduction 
from the reaction surface to the heat sink (cold rib shoulder). 
When the value of Bi/(Pr-Re-R,) increases to 6.63 (Fig. 7(b)), 
both the fluid-phase temperature and the solid-phase temper- 
ature decrease. This is because an increase in the value of 
Bi/(Pr-Re-R,) enlarges the heat transfer gate between the solid 
phase and the fluid phase that enhances the cooling of the solid 
matrices by convection. That is the hot solid matrices can trans- 
fer a large portion of heat to the cold fluid via convection. In 
addition, the entire cross section of the module middle (X = 0.8) 
has the lowest values of solid-phase temperature as the value of 
Bi/(Pr-Re-R,) increases to 6.63. 

Fig. 8(a and b) show the distributions of the solid-phase 
temperature and the fluid-phase temperature along the axial 
direction (X) at several elevations (Y=0, 0.1, and 0.3 for the 
anode, and 0.7, 0.9 and 1.0 for the cathode), respectively. The 
Re and Bi/(Pr-Re-R,;) are fixed at 6 and 6.63, respectively. As 
shown in this figure, the cathode has a higher solid-phase tem- 
perature than the anode at the corresponding location in respect 
to the module mid-plane (Y=0.5). Again, it is because of a 
significant heat generation by the electrochemical reaction in 
the cathodic catalyst layer. The concave shape for the solid- 
phase temperature distributions reflects conduction-dominant 
transport phenomena in the solid phase heat transfer. In con- 
trast, the fluid-phase temperature distribution near the shoulder 
surfaces increases downstream monotonously. At the elevation 
closed to catalyst layers, Y=0.3, and 0.7, the distribution of 
fluid-phase temperature becomes concave. Fig. 9 further shows 
the temperature distributions along the two sides of the catalyst 
layers, i.e., Y = 0.4, 0.45, 0.55 and 0.6. The temperature distribu- 


tions decrease first and then increases to a local maximum at the 
downstream end. A local minimum occurs at the module middle 
(X= 0.8). They are roughly symmetric about the local minimum 
at module middle (Y= 0.8) where faces the cold shoulder sur- 
faces. Again, the concave of the temperature distributions gives 
an indication that the conduction dominates the thermal trans- 
port in the porous medium. 

Fig. 10 shows the effect of Eco, on the temperature 
distributions at the module middle (i.e., X=0.8). The non- 
dimensional interfacial heat transfer coefficient is fixed at 
Bil(Pr-Re-R;) =6.63. As given in Table 2, an increase in Eco, 
means the increase in the overpotential heating during the elec- 
trochemical reaction. It is seen that, both the solid-phase temper- 
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Fig. 9. Temperature distributions along the two sides of catalyst layers for 
Bil(Pr-Re-R,) = 0.63. 
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Fig. 10. Effect of overpotential heating on the temperature distributions inside 
the five-layer MEA, Bi/(Pr-Re-R;,) = 6.63. 


ature and the fluid-phase temperature increase with increasing 
Eco,. Since the heat dissipation in the cathodic catalyst layer 
is more than that in the anodic catalyst layer, both of the fluid- 
phase temperature and the solid-phase temperature are higher 
in the cathodic gas diffusion layer than those in the anodic gas 
diffusion layer. 


4. Conclusions 


A numerical model has been performed to simulate the 
thermal-fluid transport phenomena in a five-layer MEA of a 
PEM fuel cell, which is in contact with interdigitated gas dis- 
tributors. Parametric studies include the flow Reynolds number 
(Re), non-dimensional heat transfer coefficient (Bi/(Re-Pr-R;)), 
and electrochemical heat parameter (Eco, ). The unique of this 
model is the first implementation of the local thermal non- 
equilibrium model in the gas diffusion layers which allows for 


a more realistic variation of the thermal-fluid in the MEA of a 
PEM fuel cell. Results show that the local maximum tempera- 
ture inside the MEA occurs at about the middle of the cathodic 
catalyst layer. In addition, an increase in the non-dimensional 
heat transfer coefficient Bi/(Pr-Re-R;) reduces the temperature 
difference between the solid matrices and the reactants in the 
gas diffusion layer. Moreover, both fluid-phase temperature and 
the solid-phase temperature decrease with increasing the non- 
dimensional heat transfer coefficient Bi/(Re-Pr-R,). It is further 
found that an increase in the overpotential heating parame- 
ter Eco, increases the fluid and solid-phase temperatures. The 
present model has successfully predicted the heat/mass transfer 
mechanisms together with the thermal pathways in the MEA 
of PEM fuel cells. It would be beneficial for further accurate 
analyses of the fuel cell thermal performance by considering the 
thermal dependent physical properties inside the MEA. 
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